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Abstract 

Modifications  of  the  Euclidean  algorithm  are  presented  for  determining  the  period  from  a  sparse 
set  of  noisy  measurements.  The  elements  of  the  set  are  the  noisy  occurrence  times  of  a  periodic 
event  with  (perhaps  very  many)  missing  measurements.  This  problem  arises  in  radar  pulse  repe¬ 
tition  interval  (PRI)  analysis,  in  bit  synchronization  in  communications,  and  other  scenarios.  The 
proposed  algorithms  are  computationally  straightforward  and  converge  quickly.  A  robust  version 
is  developed  that  is  stable  despite  the  presence  of  arbitrary  outliers.  The  Euclidean  algorithm 
approach  is  justified  by  a  theorem  which  shows  that,  for  a  set  of  randomly  chosen  positive  integers, 
the  probability  that  they  do  not  all  share  a  common  prime  factor  approaches  one  quickly  as  the 
cardinality  of  the  set  increases.  In  the  noise-free  case  this  implies  convergence  with  only  ten  data 
samples,  independent  of  the  percentage  of  missing  measurements.  In  the  case  of  noisy  data  sim¬ 
ulation  results  show,  for  example,  good  estimation  of  the  period  from  one  hundred  data  samples 
with  fifty  percent  of  the  measurements  missing  and  twenty  five  percent  of  the  data  samples  being 
arbitrary  outliers. 


1  Introduction:  Statement  of  the  Problem 


Our  problem  begins  with  a  set  of  measurements  of  a  periodic  process.  We  wish  to  determine  the 
period  of  the  process  from  this  set,  which  we  shall  assume  is  (perhaps  very)  sparse  and  noisy.  We 
assume  our  data  is  a  finite  set  of  real  numbers 

S  =  {sj}"=1 ,  with  Sj  =  kjT  +  (f>  +  rjj,  (1) 

where  r  (the  period)  is  a  fixed  positive  real  number,  the  kj's  are  non-repeating  positive  integers,  r/> 
(the  phase)  is  a  real  random  variable  uniformly  distributed  over  the  interval  [0,  r),  and  the  rjj's  are 
zero-mean  independent  identically  distributed  (iid)  error  terms.  We  assume  that  the  rj3's  have  a 
symmetric  probability  density  function  (pdf),  and  that  \rfj\  <  \  for  all  j.  By  assuming  that  the  s/s 
are  increasing,  we  can  interpret  the  s/s  as  occurrence  times  with  time  gaps  or  jumps  determined 
by  the  kj's.  (Later  we  will  have  occasion  to  reorder  the  data  for  analysis;  this  does  not  alter  the 
time-of-occurrence  interpretation).  For  example,  letting  k\  =  1,  =  3,  £3  =  5,  . . .,  kn  =  2n  —  1, 

corresponds  to  missing  every  other  occurrence  of  a  periodic  event  with  period  r.  The  period  is  the 
unique  maximum  value  of  r  such  that  the  kj  are  all  integers.  Note  that  for  any  solution  t  that 
fits  the  form  of  S,  the  numbers  |,  |, . . .  are  also  possible  solutions.  Estimates  of  r  in  turn  lead  to 
estimates  of  (j)  and  the  kj's. 

The  problem  of  finding  r  given  the  set  (1)  arises  in  different  ways  that  are  often  associated 
with  pulse  train  analysis.  These  include  radar  pulse  repetition  interval  (PRI)  analysis  [31],  bit 
synchronization  in  communication  systems  [10],  and  neuron  firing  rate  analysis.  More  generally, 
the  problem  arises  when  one  observes  events  that  occur  in  integer  multiples  of  some  fixed,  but 
unknown,  real  number.  The  problem  is  to  determine  that  number.  For  example,  in  communications 
the  elements  of  S  may  be  the  times  of  a  pseudorandomly  occurring  change  in  the  carrier  frequency, 
where  the  change  rate  is  governed  by  a  shift  register  output.  The  set  S  may  also  arise  from  noisy 
zero-crossing  times  of  a  sinusoid. 

We  solve  for  r  by  using  a  modified  Euclidean  algorithm  to  find  what  is  essentially  the  greatest 
common  divisor  of  the  set  S.  Subsequent  algorithms  are  modifications  of  this  first  procedure.  The 
theoretical  basis  for  the  algorithm  is  given  in  Sections  2,  3,  4,  and  the  Appendix.  The  algorithm 
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is  computationally  straightforward,  requires  little  input  data,  and  given  reasonable  data,  quickly 
converges  to  the  value  of  r  with  very  high  probability.  A  fundamental  result  is  given  by  Theorem  3.1, 
which  states  that  for  a  set  of  uniformly  distributed  positive  integers  distributed  over  an  arbitrarily 
large  lattice,  the  probability  that  they  all  do  not  share  a  common  prime  factor  approaches  one 
quickly  as  the  cardinality  of  the  set  increases. 

In  the  noise-free  case  our  algorithm  is  equivalent  to  the  Euclidean  algorithm  and  converges 
with  very  high  probability  given  only  n  =  10  data  samples,  independent  of  the  number  of  missing 
measurements.  Simulation  examples  demonstrate  successful  estimation  of  r  for  n  =  10  with  99.99% 
of  the  possible  measurements  missing.  In  fact,  with  only  10  data  samples,  it  is  possible  to  have 
the  percentage  of  missing  measurements  arbitrarily  close  to  100%.  There  is,  of  course,  a  cost,  in 
that  the  number  of  iterations  the  algorithm  needs  to  converge  increases  with  the  percentage  of 
missing  measurements.  The  algorithm  will  degrade  as  n  is  reduced  from  10  (see  Table  2),  which  is 
consistent  with  the  theory  (Theorem  3.1,  Proposition  3.1,  and  Table  1). 

In  the  presence  of  noise  (non-zero  rjj's  in  (1))  and  false  data  (or  outliers),  there  is  a  tradeoff 
between  the  number  of  data  samples,  the  amount  of  noise,  and  the  percentage  of  outliers.  The 
algorithm  will  perform  well  given  low  noise  for  n  =  10,  but  will  degrade  as  noise  is  increased. 
However,  given  more  data,  it  is  possible  to  reduce  noise  effects  and  speed  up  convergence  by 
binning  the  data,  and  averaging  across  bins.  Binning  can  be  effectively  implemented  by  using  an 
adaptive  threshhold  with  a  gradient  operator,  allowing  convergence  in  a  single  iteration  in  many 
cases. 

The  problem  of  finding  r  from  S  has  been  considered  by  Fogel  and  Gavish  [9],  who  used 
spectrum  analysis  of  a  delta-train  function  with  impulses  at  the  occurrence  times.  This  approach 
is  equivalent  to  a  periodogram-based  method  for  spectral  analysis  of  point  processes  as  developed 
by  Bartlett  [1].  General  techniques  for  PRI  analysis  include  dividing  the  total  number  of  pulses 
observed  by  the  number  of  pulse  intervals,  as  well  as  use  of  pulse  interval  histograms  [31].  More 
recently  Gray  et  al.  have  derived  the  maximum  likelihood  estimator  of  t  for  the  iid  Gaussian  noise 
case  [12].  In  their  work,  missing  measurements  and  outliers  were  not  considered. 

The  more  general  problem  of  time  series  spectrum  analysis  using  zero-crossings  has  been  con- 
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sidered  by  a  number  of  authors.  Kedem  uses  zero-crossing  counts  on  an  interval  for  spectrum 
analysis  of  random  signals  [18].  This  method  does  not  seem  to  be  conducive  to  cases  of  missing 
measurements.  Jones  [16]  and  Parzen  [22]  developed  a  multiplicative  missing  observation  model 
(AM  model)  for  spectrum  analysis  of  stationary  processes  with  periodically  missing  measurements. 
The  resulting  nonparametric  estimates  are  based  on  suitably  modified  periodograms.  This  work 
was  extended  by  Bloomfield  [2]  and  Scheinok  [29]  to  include  missing  measurements  where  the  prob¬ 
ability  of  a  miss  is  modeled  by  a  Bernoulli  process.  Identification  of  parametric  (ARMA)  models 
from  time  series  with  missing  measurements  has  also  been  considered  by  a  number  of  authors, 
e.g.,  see  Rosen  and  Porat  [27]  and  references  therein.  Recently  Dandawate  [6]  and  Giannakis  and 
Zhou  [11]  have  proposed  using  cyclic  statistics  of  time  series  with  periodically  missing  measure¬ 
ments.  Generally  speaking,  the  above  referenced  models  for  missing  measurements  assume  access 
to  uniformly  spaced  amplitude  samples,  and  are  therefore  not  directly  applicable  in  our  setup. 

Kay  and  Sudhaker  proposed  the  use  of  the  occurrence  time  of  zero-crossings  to  obtain  DFT 
coefficients  [17].  This  approach  is  related  to  a  class  of  pulse  time  modulation  techniques  in  which 
sample  amplitudes  are  suitably  transformed  to  zero-crossing  interval  widths  [8].  This  method 
requires  addition  of  an  auxiliary  signal  before  detection  of  the  zero-crossings  and  is  therefore  not 
applicable  here. 

The  paper  is  organized  as  follows.  Section  2  gives  the  development  of  a  modified  Euclidean 
algorithm  for  finding  r  that  is  robust  to  noise.  Given  this  estimate,  we  then  discuss  estimating 
<f)  and  the  kj's.  Section  3  is  a  discussion  of  the  result  that  for  a  set  of  randomly  chosen  positive 
integers,  the  probability  that  they  do  not  all  share  a  common  prime  factor  approaches  one  quickly 
as  the  cardinality  of  the  set  increases.  We  then  propose,  in  Section  4,  further  modifications  of 
the  algorithm  which  are  effective  given  increased  noise  and  outliers.  Section  5  contains  simulation 
results  and  discussion  of  performance  for  the  various  proposed  algorithms.  Lengthy  proofs  are 
presented  in  the  Appendix. 


3 


2  A  Modified  Euclidean  Algorithm  for  Finding  r 


Given  the  set  S  as  in  (1),  we  modify  the  Euclidean  algorithm  to  develop  a  procedure  for  finding 
t.  This  involves  some  basic  number  theory.  References  for  this  material  include  Hardy  and  Wright 
[13],  Ireland  and  Rosen  [14],  Knuth  [19],  Leveque  [21],  Rosen  [26],  and  Schroeder  [30]. 

The  Euclidean  algorithm  is  a  division  process  for  the  set  of  integers  Z.  The  algorithm  is  based 
on  the  property  that,  given  two  positive  integers  a  and  b,  a  >  b,  there  exists  two  positive  integers 
q  and  r  such  that 

a  =  q-  b  +  r,0<r<b. 

If  r  =  0,  we  say  that  b  divides  a,  and  denote  this  by  b\a.  This  property  of  the  set  of  integers, 
combined  with  the  fact  that  if  a,  6  €  Z  \  {0}  then  a  ■  b  ^  0  (Z  has  no  zero  divisors),  make  Z  a  unique 
factorization  domain.  Thus  in  Z  every  non-zero  element  may  be  written  as  the  product  of  powers 
of  irreducible  integers,  or  primes.  The  Euclidean  algorithm  also  holds  in  more  general  algebraic 
structures  called  Euclidean  domains. 

The  Euclidean  algorithm  yields  the  greatest  common  divisor  of  two  (or  more)  elements  of  Z. 
The  greatest  common  divisor  of  two  integers  a  and  b,  denoted  by  gcd(a,  b),  is  the  product  of  the 
powers  of  all  prime  factors  p  that  divide  both  a  and  b.  We  may  represent  the  algorithm  applied  to 
a,b,  a>  b,  as  follows: 

a  =  b  ■  q\  +  ri  :  0  <  r\  <  b 

b  =  r\  •  <72  +  r2  :  0  <  r2  <  r\ 

n-2  =  rk- 1  -qk  +  rk  ■  0  <  rk  <  rk-i 

n- i  —  T'k'  qk- 


The  procedure  terminates  when  rk+ 1  =  0.  This  gives  gcd(a,  b)  =  rk- 

This  procedure  can  be  extended  to  work  on  S.  The  symbol  gcd(Aq, . . . ,  kn)  is  the  greatest 
common  divisor  of  the  set  {kj},  i.e. ,  the  product  of  the  powers  of  all  prime  factors  p  that  divide 
each  kj.  Note  that  this  is  not  the  pairwise  gcd  of  the  set  {kj}.  If  gcd(Aq, . . . ,  kn)  =  1,  the  set 
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{kj}  is  called  mutually  relatively  prime.  If,  however,  gcd(ki,kj)  —  1  for  all  i  /  j,  the  set  {kj}  is 
called  pairwise  relatively  prime.  If  a  set  is  pairwise  relatively  prime,  it  is  mutually  relatively  prime. 
However,  the  converse  is  not  true  (for  example,  consider  the  set  {35,21, 15}).  The  computation  of 
the  gcd  of  a  set  of  more  than  two  integers  uses  the  following  proposition.  There  is  also  a  natural 
extension  of  the  gcd  to  multiples  of  a  fixed  r  >  0. 

Proposition  2.1 

(*■)  gcd(&ir, . . . ,  knr)  =  Tgcd(fci,...,fcn), 

(**■)  gcd  (ki,  ...,kn)  =  gcd(Aq, . . . ,  fcn_2,  (gcd(A:n_1,  kn))) . 

Proof  :  See  Leveque  [21],  page  16.  □ 

The  standard  Euclidean  algorithm,  as  shown  above,  involves  repeated  division.  In  our  problem, 
we  are  dealing  with  numbers  that  are  essentially  “noisy  integers.”  Remainder  terms  could  be 
noise,  and  thus  could  be  non-zero  numbers  arbitrarily  close  to  zero.  Subsequent  iterations  in 
the  procedure  may  involve  dividing  by  such  small  values,  which  would  result  in  arbitrarily  large 
numbers.  Thus,  the  standard  algorithm  is  unstable  under  perturbation  by  noise.  However,  the 
algorithm  may  be  changed  so  that  the  process  of  subtraction  replaces  division  by  making  use  of  the 
following  proposition.  So  that  (kj  —  kj+i)  €  N,  we  assume  that  the  kj: s  are  sorted  in  descending 
order. 

Proposition  2.2 

gcd(fci, . . . ,  kn)  =  gcd((Aq  -k2),(k2  -  fc3), . . . ,  (fc„_ i  -kn),kn). 


Proof  :  See  the  Appendix.  □ 

We  will  also  need  the  following  related  result. 

Proposition  2.3 

gcd((fci  -  k2),  (k2  -  k3), ...,  (kn-i  -  kn))  =  gcd((A:1  -  kn),  (k2  -  kn), ...,  (kn- 1  -  kn)) . 
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Proof  :  See  the  Appendix.  □ 


We  are  now  ready  to  develop  our  modified  algorithm.  Here,  and  throughout  the  rest  of  the 
paper,  we  will  assume  the  s/s  are  sorted  in  descending  order,  i.e.,  S\  >  S2  >  ■  •  •  >  sn.  This  allows 
a  more  straightforward  visualization  of  our  algorithm.  We  form  a  new  set  by  subtracting  adjacent 
pairs  of  these  numbers,  given  by  Sj  —  Sj+i.  After  this  first  operation,  the  phase  information  has 
been  subtracted  out,  and  the  resulting  set  has  the  simpler  form 

S'  =  ,  with  s'  =  K3t  +  77'  , 

where  Kj  =  kj  —  kj+\  and  77'  =  777  —  777+1.  In  subsequent  iterations  of  the  algorithm,  the  data 
will  maintain  this  same  general  form.  Because  of  the  777  perturbations  we  establish  a  threshold 
770  and,  after  the  subtraction  of  adjacent  pairs,  we  declare  all  numbers  in  the  interval  [0, 770]  to  be 
zero  and  eliminate  them  from  the  set.  Choice  of  770  is  dictated  by  the  distribution  of  the  777*8,  with 
0  <  770  <  §  (see  Section  4).  We  then  append  zero,  sort,  subtract  adjacent  pairs,  and  threshold.  By 
appending  zero,  we  adjoin  the  previous  non-zero  minimum  to  the  set.  The  algorithm  is  continued 
by  iterating  this  process  of  appending  zero,  sorting,  subtracting,  and  eliminating  the  elements  in 
[0,770].  It  terminates  when  all  but  one  of  the  elements  are  in  [0,770],  i.e., “equal  to  zero.”  By  the 
propositions  above,  this  element  is  equal  to  gcd(Ai, . . . ,  i£n-l)  *  T  ±  error  term . 

We  will  see  in  the  next  section  that  gcd(A'i, . . . ,  ATn-i)  — >  1  with  probability  1  as  n  — >  00. 
Moreover,  we  will  see  that  this  convergence  is  very  fast,  which  shows  that  the  proposed  modified 
Euclidean  algorithm  yields  r  with  a  high  probability  for  small  (n  ~  10)  to  moderate  (n  ~  100) 
values  of  n,  depending  on  the  distribution  of  the  777 ’s  and  the  percentage  of  outliers. 

The  algorithm  for  finding  r  given  the  set  S  is  summarized  as  follows.  Again,  we  assume  that 
the  original  data  set  is  in  descending  order,  i.e.,  Sj  >  Sj+ 1. 

Modified  Euclidean  Algorithm 


1. )  After  the  first  iteration,  append  zero. 

2. )  Form  the  new  set  with  elements  Sj  —  S7+1. 
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3. )  Sort  in  descending  order. 

4. )  Eliminate  elements  in  [0,  r/o]  from  end  of  the  set. 

5. )  Algorithm  is  done  if  left  with  a  single  element.  Declare  t  =  sj.  If  not  done,  go  to  (1.). 


Remarks  : 


1.  It  is  easily  possible  to  generate  low  noise  or  even  noise-free  data  sets  such  that  the  modified 
Euclidean  algorithm  does  not  converge  to  r,  even  for  an  arbitrarily  large  number  of  data 
points.  For  example,  consider  k\  =  3,  —  5,  —  7, . . .  with  true  underlying  period  r  =  1. 

To  eliminate  the  phase,  we  do  not  save  the  first  minimum.  For  this  data,  the  algorithm  will 
give  t  —  2.  However,  for  even  10  data  points,  the  theory  developed  in  the  following  section 
shows  that  such  data  sets  are  extremely  unlikely  to  occur.  Another  way  of  saying  this  is  that 
if  a  large  set  of  low  noise  or  noise-free  data  leads  the  algorithm  to  an  incorrect  solution,  then 
the  corresponding  k^s  exhibit  a  consistent  pattern. 


2.  Given  an  estimate  f  of  r,  it  is  then  possible  to  estimate  p  and  the  kj’s.  We  give  a  brief  dis¬ 
cussion  of  these  estimates.  Fogel  and  Gavish  [9,  10]  have  addressed  the  problem  of  estimating 
(f>,  and  have  produced  an  estimator  p  for  Pi  given  by 


y~]  exp(27n 


J= 1 


(2) 


For  a  good  estimate  f  of  r,  exp(2ni^)  ~  1.  For  rjj  <C  r  «  r,  r/j/f  1,  and  so 
exp(27ru/?/f)  «  exp(0)  =  1.  Therefore, 


^  I  n  s 

—  arg  <  ^exp(27ri^) 


LJ=1 

n 


=  ~  arg  |y^exp(27ri^-)exp(27r^)exp(2-7r^) 

T  \  n  (f) 

EexP(2«|) 


j=i 


=  arg  in  •  exp(27rj^)l 
Z7T  t  T  J 
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r 

27T 


arg  {n}  +  arg  |exp(27ri^)|^ 
=  £arg{exp(2™|)} 


T  27 T(j) 
2i r  f 


=  <t> 


3.  We  present  two  methods  of  getting  an  estimate  on  the  set  of  kj’s.  Given  a  good  estimate  <fi, 
the  first  is  to  form  the  set 

a  =  {kjT  +  <j>  +  7]j-  4>}nj=l . 

Given  the  estimate  f,  estimate  kj  by 

%  =  round  (kL+l±J!Lz£\  ,  (3) 

where  round(-)  denotes  rounding  to  the  nearest  integer. 

We  could  also  work  with  the  set 


o'  =  iKiT  +  Vj}]=i  U  {knr  +  <p  +  r]n-$}, 
where  Kj  =  kj  —  k]+\  and  ry'  =  rjj  —  r)j+ Given  the  estimate  f,  estimate  kn  by 

kn  =  round 


and  Kj  by 


Kj  =  round  ( 


KjT  +  rj'j 


)• 


Then,  kn-i  =  Kn- 1  +  kn,  kn-2  =  Kn- 2  +  fcn-i,  and  so  on. 


3  Relatively  Prime  Numbers  and  the  Zeta  Function 

In  this  section,  we  present  Theorem  3.1,  which  says  that  gcd(&i, ...  ,kn)  — >  1  with  probability  1 
as  n  — >  00.  Moreover,  this  convergence  is  very  quick.  We  use  this  to  show  that 

gcd(KTiT,  •  •  • ,  Kn-\T) — >  r,  (4) 
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with  probability  1  as  n  — >  oo.  We  show  that  the  proposed  modified  Euclidean  algorithm  yields 
r  with  a  high  probability  for  small  (n  «  10)  to  moderate  (n  fts  100)  values  of  n,  depending  on  the 
distribution  of  the  ly’ s  and  the  percentage  of  outliers.  In  the  noise-free  case,  the  theory  tells  us 
that  the  algorithm  will  almost  certainly  yield  r  given  only  10  data  samples. 

The  proof  is  relatively  long,  and  is  presented  in  the  Appendix.  The  development  of  the  theorem 
involves  techniques  and  results  from  analytic  number  theory.  Two  functions,  Riemann’s  Zeta 
function  ((n)  and  the  Mobius  function  /i(m),  play  key  roles.  The  result  is  classical  for  the  case 
n  —  2,  and  was  proven  by  Dirichlet  in  1849  (see  Knuth  [19],  pp.  324,  337,  595,  and  Schroeder  [30], 
pp.  48-50).  An  elegant  formulation  for  n  =  2  is  given  in  Hardy  and  Wright  [13]. 

Riemann’s  Zeta  function  is  defined  on  the  complex  half  space  {z  G  C  :  !?(z)  >  1}  by 

OO 

C(z)  =  Y,n~z-  (5) 

71—  1 

Euler  demonstrated  the  connection  of  £  with  number  theory  by  showing,  in  1736,  that 

oo  , 

where  P  =  {pi,P2,Pz,  •  •  ■}  =  {2, 3, 5, . . .}  is  the  set  of  all  prime  numbers  (see  Conway  [5]).  In  the 
following  P{-}  denotes  probability. 

Theorem  3.1  Given  n  (n  >  2)  randomly  chosen  positive  integers  {k\, . . .  ,kn}, 

P{gcd(fci,...,fc„)  =  1}  =  [C(n)]_1  • 

The  proof  of  Theorem  3.1  follows  directly  from  the  following  theorem.  We  let  card{-}  de¬ 
note  cardinality  of  the  set  {•},  and  let  {1, . . .  ,£}n  denote  the  sublattice  of  positive  integers  in  Rn 
with  coordinates  c  such  that  1  <  c  <  l.  Therefore,  Nn(t)  =  card{(Aq, . . . ,  kn)  e  {l,...,£}n  : 
gcd(Aq, . . . ,  kn)  =  1}  is  the  number  of  relatively  prime  elements  in  {1, . . . 

Theorem  3.2  Let 


Nn(i)  =  card{(fci,...,fc„)  e  {1  ,...,£}n  :  gcd{ku...  ,kn)  =  1}  , 


For  n  >  2,  we  have  that 


lim 

t — >oo 


Nn(£) 

Hjn 


(C(n)]"1  - 
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Proof  :  See  the  Appendix.  □ 


Our  original  data  set  contains  phase  information.  We  eliminate  this  by  the  subtraction  of 
adjacent  elements  in  the  set.  But  then,  rather  than  working  with  {&i, . . . ,  kn},  we  are  working  with 
{(&!  —  /c'2) , . . . ,  (kn- 1  —  kn)}.  Therefore,  we  need  the  following. 

Corollary  3.1  Let  {k\, . . .  ,kn}  be  n  (n  >  3)  randomly  chosen  positive  integers,  with  kj  >  kj+ 
and  let  Kj  =  kj  —  kj+i  for  j  =  1, . . . ,  n  —  1.  Then 

P{gcd(Ku. . . ,  Kn-!)  =  1}  =  [C(n  -  I)]’1  . 


Proof  :  See  the  Appendix.  □ 


Various  techniques,  e.g.,  Cauchy  residue  theory,  allow  us  to  get  a  closed  form  solution  of  ((n) 
for  the  even  natural  numbers.  This  is  (see  Ireland  and  Rosen  [14]  and  Conway  [5]) 

\2  k 


2C(2A)  =  Bernoulli  Number)  . 

(2k)\ 


(7) 


(“The  result  is  due  to  Euler  and  constitutes  one  of  his  most  remarkable  computations.”  Ireland  and 
Rosen  [14],  page  231.)  The  values  ((2k  +  1)  can  be  estimated  numerically  (see  the  CRC  tables  [3]). 
Some  representative  values  are  shown  in  Table  1.  We  can  see  from  these  values  that  [((n)]”1  — ->  1 
quickly  as  n  increases.  In  fact,  we  can  show  that  this  rate  of  convergence  is  at  least  exponential. 


Proposition  3.1  Let  co  €  (l,oo).  Then 

Jim  [CMP1  =  i, 

converging  to  1  from  below  faster  than  1/(1  —  21-w). 


(8) 


Proof  :  See  the  Appendix.  □ 

Combining  Corollary  3.1  with  Propositions  2.1  -  3.1  gives  the  theoretical  underpinnings  of  the 
algorithm. 
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Theorem  3.3  Let  {ki,...,kn}  be  n  (n  >  3)  randomly  chosen  positive  integers,  with  kj  >  kj+\, 
and  let  Kj  =  kj  —  kj+\  for  j  =  1, . . .  ,  n  —  1.  Then 

gcd(ifir, . . . ,  Kn-ir)  — »  r , 

with  probability  1  as  n  — >■  oo. 


We  close  this  section  with  a  heuristic  argument  for  the  validity  of  Theorem  3.1.  The  argument, 
although  not  rigorous,  does  provide  insight  into  the  result.  In  this  discussion,  we  have  not  carefully 
defined  what  it  means  for  the  integers  to  be  chosen  at  random,  and  we  do  not  present  the  estimates 
needed  to  justify  convergence.  These  are  presented  in  the  proof  of  Theorem  3.2  in  the  Appendix. 

Given  randomly  distributed  positive  integers,  by  the  Law  of  Large  Numbers,  about  1/2  of  them 
are  even,  1/3  of  them  are  multiples  of  three,  and  1/p  are  a  multiple  of  some  prime  p.  Thus,  given 
n  independently  chosen  positive  integers, 

P{p\ki,p\k2, ,  and  p\kn}  = 

P{p\h}  ■  P{p\k2}  ■  P{p\kn}  = 

1/(P)  ■!/(?)••■  ■■!/(?)  = 

1  /(P)n  • 


Therefore, 


P{P  Xh,P  J(k2,  ■  ■  •  ,or  p  J(kn}  =  1  -  l/(p)n  . 


(9) 


By  the  Fundamental  Theorem  of  Arithmetic,  every  integer  has  a  unique  representation  as  a  product 
of  primes.  Combining  that  theorem  with  the  definition  of  gcd,  we  get 

OO 

P{gcd(A;i, . . . ,  kn)  =  1}  =  IJ  1  -  1  /{pj)n.  (10) 

3= 1 

where  pj  is  the  jth  prime.  But,  by  Euler’s  formula, 

OO  -t 

cw-n-— (ii) 

Therefore, 


P{gcd{ki,...  ,kn)  =  1}  =  l/(C(n)). 


(12) 


This  completes  our  heuristic  argument. 
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4  Further  Modifications 


In  this  section  we  consider  the  effects  of  noise  and  outliers  on  the  modified  Euclidean  algorithm 
described  in  Section  2.  The  modified  algorithm  replaces  division,  as  required  in  the  standard 
Euclidean  algorithm,  with  repeated  subtraction  in  order  to  gain  stability  with  respect  to  noise. 
Error  analysis  of  this  approach  is  complicated  by  the  facts  that  the  algorithm  is  iterative  and  that 
the  algorithm  sorts  the  data.  Therefore,  this  analysis  involves  order  statistics. 

Suppose  the  pdf  of  the  r/j's  is  given  by  /^(y),  and  consider  the  set  of  differences  obtained  in  the 
first  iteration,  given  by 

Vj  =  Sj  -  8j+ 1  =  {kj  -  kj+\)r  +  (rjj  -  rjj+x).  (13) 

Invoking  the  zero-mean  iid  assumption  on  the  rjj's,  the  pdf  of  ( rlj~rlj+i )  is  given  by  the  convolution 
frjiv)  *  ftjiv)-  So,  for  example,  if  fv{rj)  ~  U[—  y,  y]  (t?  is  uniformly  distributed  with  parameter  A) 
then  fy  (y)  =  tri[y—  (kj — kj+i)r],  the  triangle  function  centered  at  (kj—kj+i)r.  Two  points  can  now 
be  made.  First,  after  the  first  iteration,  the  differencing  operation  has  removed  the  independence  of 
the  error  terms.  Second,  the  ordering  operation  makes  the  nature  of  the  dependence  in  subsequent 
iterations  difficult  to  determine.  Analysis  of  order  statistics  very  often  rests  on  an  iid  assumption, 
e.g.,  see  Sarhan  and  Greenberg  [28]  and  Reiss  [25].  Without  the  iid  assumption,  this  analysis  leads 
into  many  open  questions  (see  Reiss  [25]). 

In  general,  beyond  the  first  iteration  the  pdf  of  the  subsequent  error  terms  becomes  asymmetric, 
even  when  starting  with  iid  rjj' s  with  symmetric  pdf  fTI  (rj) .  This  occurs  due  to  the  reordering  before 
differencing  at  each  iteration,  and  because  after  the  first  iteration  the  errors  are  no  longer  iid.  The 
result  is  that  using  the  modified  Euclidean  algorithm  can  lead  to  negatively  biased  estimates  of  r 
after  the  first  iteration  due  to  the  skewness  of  the  pdf  of  the  errors.  As  we  will  see  this  can  be 
corrected  for  by  averaging. 

In  order  to  illustrate  the  behavior  of  the  algorithm  consider  the  following  example.  Let  the  set 
S  of  equation  (1)  be  generated  as  follows.  Let  r  =  1,  n  =  100  data  samples,  the  jumps  in  the  kj's 
be  randomly  selected  from  a  discrete  uniform  distribution  on  the  interval  [1, 10],  and  the  noise  be 
iid  and  uniformly  distributed  as  fv{rj)  ~  U[— 0.1, 0.1].  A  data  set  S  was  generated  according  to 
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these  parameters  and  used  as  input  to  our  algorithm.  Consider  the  results  after  one  iteration,  in 
which  the  data  has  been  differenced  and  sorted  into  descending  order,  as  plotted  in  Figure  1.  The 
data  are  clustered  into  “steps”  around  integer  multiples  of  r  =  1,  as  we  expect  from  (13).  That  the 
steps  are  all  of  the  same  approximate  length  is  due  to  the  uniform  distribution  in  the  jumps  of  the 
kj's  in  the  original  data  set  5.  Other  distributions  will  result  in  different  proportions.  From  (13) 
and  the  assumptions  on  the  noise  we  know  that  the  data  has  a  mean  that  is  an  integer  multiple 
of  r  given  by  ( kj  —  kj+i)r,  with  noise  symmetrically  distributed  around  this  mean.  This  suggests 
isolating  each  step  and  averaging  the  data  within  each  step  to  reduce  noise  effects. 

A  straightforward  method  for  clustering  the  data  is  to  employ  a  gradient  operator  to  determine 
when  a  step  has  occurred.  After  the  first  iteration  (as  in  Figure  1)  the  gradient  is  estimated,  with 
large  gradient  values  indicating  a  step  or  “edge”  in  the  data.  We  have  employed  a  simple  estimator 
by  convolving  with  an  impulse  response  given  by  [-1,0, 1].  This  operator  is  well  known  in  signal 
and  image  processing,  e.g.,  see  Jain  [15].  A  data-adaptive  gradient  threshold  go  is  selected  as  10% 
of  the  maximum  gradient  value,  and  data  points  above  this  threshold  are  assumed  to  correspond 
to  the  step  edges.  After  the  steps  have  been  isolated  the  step  heights  are  easily  found,  and  the 
minimum  step  height,  call  it  f,  is  taken  as  a  coarse  estimate  of  r.  Referring  to  Figure  1,  all  of 
the  step  heights  are  approximately  equal  to  r,  again  due  to  the  original  distribution  of  the  jumps 
in  the  kj' s  used  in  generating  S.  We  then  use  f  to  set  two  thresholds.  The  first  is  rjo  =  0.35f, 
used  to  define  the  neighborhood  of  zero  in  which  data  will  be  eliminated  during  each  iteration  (see 
Section  2).  The  second  we  take  to  be  yo  =  0.6f,  used  to  segment  the  steps  at  each  iteration.  The 
segmentation  proceeds  by  searching  for  jumps  in  height  greater  than  yo,  and  averaging  over  each 
segment.  The  choices  of  0.35  and  0.6  are  based  on  extensive  simulation  experience  (see  Section  5), 
and  can  be  more  rigorously  justified  in  specific  cases.  However,  performance  is  reasonably  robust  to 
changes  in  these  weights  under  the  various  scenarios  considered.  The  averaging  produces  significant 
data  reduction,  and  therefore  greatly  increases  the  speed  of  convergence.  The  gradient  operator  is 
applied  only  as  part  of  the  first  iteration,  the  data  reduces  rapidly  with  each  iteration  and  precludes 
use  of  the  gradient  operator  except  as  part  of  the  first  iteration. 

We  summarize  the  foregoing  in  the  following  algorithm  statement.  Again  we  assume  the  data 
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is  initially  sorted  in  descending  order.  Recall  that  appending  zero  in  the  first  step  appends  the 
previous  minimum. 


Modified  Euclidean  Algorithm  (With  Averaging) 


1. )  After  the  first  iteration,  append  zero. 

2. )  Form  the  new  set  with  elements  Sj  —  Sj+i. 

3. )  Sort  in  descending  order. 

4. )  On  the  first  iteration,  apply  gradient  and  obtain  f,  yielding  r?o  =  0.35t  and  yo  =  0.6f  (see 

text). 

5. )  Average  the  data  over  each  step,  with  steps  determined  by  jumps  of  height  yo- 

6. )  Eliminate  elements  in  [0,  r/o]  from  end  of  the  set. 

7. )  Algorithm  is  done  if  left  with  a  single  element.  Declare  t  =  s\.  If  not  done,  go  to  (1.). 

If  the  data  is  such  that,  after  the  first  iteration,  there  is  a  relatively  large  cluster  around  the  step 
nearest  to  zero,  then  we  can  readily  estimate  r  by  finding  this  step,  averaging  only  over  these  data 
points,  and  declaring  this  to  be  t.  This  is  the  rightmost  or  lowest  step  after  the  first  iteration  (see 
Figure  1).  Under  our  assumptions  this  mean  is  an  unbiased  estimate  of  r.  Accurate  estimation  of  r 
from  a  single  iteration  assumes  that  n  is  large  enough,  and  the  span  of  the  original  data  set  S  small 
enough,  to  yield  sufficient  data  in  the  neighborhood  of  t.  This  is  a  function  of  the  distribution  of 
the  kj' s. 

In  our  example  above  the  missing  observations  were  modeled  by  taking  the  jumps  in  the  kj' s 
as  uniformly  distributed  on  the  discrete  interval  [1  ,  M],  with  M  —  10.  Thus,  for  large  n,  after  the 
first  iteration  the  data  will  cluster  in  M  steps  with  an  expected  value  of  n/M  samples  in  each  step, 
as  in  Figure  1.  As  another  model  for  missing  observations  we  can  employ  an  iid  Bernoulli  process 
to  determine  if  an  observation  is  missing  or  not,  where 

P(missing  observation)  =  A 

P(observation  occurring)  =  1  —  A.  (14) 
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For  example,  with  A  =  0.6,  we  expect  60%  of  the  observations  to  be  missing.  Given  an  element  of 
S,  Sj  —  kjT  +  rjj  +  (j),  then  the  probability  that  Sj+i  =  kj+ir  +  rjj+ 1  +  cf>  =  (kj  +  l)r  +  +  (f) 

is  an  element  of  the  set  S  is  given  by  1  —  A.  It  follows  that,  after  the  first  iteration,  we  expect 
(1  —  A)(n  —  1)  data  samples  to  be  in  the  lowest  step  clustered  around  the  true  value  of  r.  Thus,  if 
n  =  101  and  A  =  0.6,  then  after  the  first  iteration  we  expect  (1  —  0.6)100  =  40  data  samples  to  be 
clustered  in  the  lowest  step  around  r;  the  average  over  these  40  data  samples  may  then  be  taken 
as  an  estimate  of  r. 

The  single  iteration  algorithm  is  a  simple  modification  of  the  preceding  multi-iteration  version. 
After  the  (first)  difference  and  sort  operations  the  gradient  is  applied  and  the  lowest  step  isolated. 
The  average  over  this  step  is  then  taken  as  our  estimate  of  r.  We  summarize  the  single  iteration 
form  of  the  algorithm  as  follows. 

Modified  Euclidean  Algorithm  (Single  Iteration) 


1. )  Given  the  set  S,  form  the  new  set  with  elements  sj  —  sJ+1. 

2. )  Sort  the  new  set  in  descending  order. 

3. )  Apply  gradient  operator  and  obtain  f,  with  yo  =  0.6f. 

4. )  Obtain  t  by  averaging  the  data  over  the  lowest  step,  isolating  this  step  based  on  the  lowest 

two  jumps  of  height  yo- 

Next  we  consider  the  effect  of  arbitrary  outliers.  These  are,  in  general,  quite  harmful  to  the 
estimation  of  the  gcd.  This  is  easily  seen  in  the  noise-free  case  because  the  gcd  of  a  contaminated 
set  may  be  arbitrarily  different  from  the  gcd  of  the  uncontaminated  set.  Testing  of  the  algorithms 
presented  thus  far  shows  sensitivity  to  the  presence  of  even  a  single  outlier.  This  is  because  the 
outliers  will  not  necessarily  fall  into  the  step-like  clusters  we  expect. 

A  more  robust  version  of  the  algorithm  can  be  obtained  with  the  introduction  of  a  step-width 
threshold,  xo.  Let  us  concentrate  on  the  single  iteration  approach.  After  application  of  the  gradient 
operator  as  before,  we  avoid  false  edges  and  clusters  occurring  below  the  true  lowest  step  by 
requiring  the  step-width  (the  number  of  data  samples  in  a  particular  step)  to  be  greater  than  xq  . 
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The  choice  of  xo  is  dictated  by  two  considerations,  the  distribution  of  the  outliers  and  the 
expected  number  of  data  samples  in  the  lowest  step  after  the  first  iteration.  This  represents  a 
tradeoff.  For  example,  consider  again  the  Bernoulli  model  for  missing  observations.  Suppose 
A  =  0.25  and  n  =  101.  We  therefore  expect  75  data  samples  in  the  lowest  step,  allowing  us  to 
choose  xq  —  50,  say.  The  single  iteration  algorithm  employing  xq  will  therefore  search  for  the  lowest 
step  in  the  data  that  contains  at  least  xo  =  50  data  samples.  Depending  on  the  distribution  and 
number  of  outliers  it  may  be  extremely  unlikely  that  a  step  will  occur  that  simultaneously  has  a 
mean  less  than  r  and  has  more  than  xq  =  50  data  values.  The  drawback  to  employing  x$  is  that 
setting  its  value  requires  some  a  priori  knowledge  or  guesswork  and  is  not  easily  set  adaptively.  The 
advantage  is  that  for  reasonable  scenarios  its  use  makes  the  algorithm  very  robust  to  the  presence 
of  even  large  percentages  of  outliers,  as  illustrated  in  the  simulation  results  of  the  next  section. 
This  robust  single  iteration  algorithm  is  summarized  as  follows. 

Robust  Single  Iteration  Algorithm 

1. )  Given  the  set  S,  form  the  new  set  with  elements  Sj  —  Sj+\. 

2. )  Sort  the  new  set  in  descending  order. 

3. )  Apply  gradient  operator  and  obtain  f,  with  yo  —  0.6f. 

4. )  Obtain  f  by  averaging  the  data  over  the  lowest  step,  isolating  the  step  based  on  the  lowest 

two  jumps  of  height  yo  with  step-width  greater  than  xo- 

Further  discussion  of  the  various  algorithms  and  their  performance  is  presented  at  the  end  of  the 
next  section. 

Remarks  : 

1.  The  computational  load  of  the  robust  single  iteration  algorithm  can  be  approximated  as 
follows.  The  differencing  operation  of  1.)  requires  n  —  1  additions,  a  fast  sort  can  be  accom¬ 
plished  in  about  nlogn  operations  [20],  and  convolution  with  the  gradient  operator  [—1,0, 1] 
can  be  accomplished  with  n  additions  and  requires  no  multiplies.  Finally,  4.)  requires  a  few 
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differences  and  comparisons.  Thus,  the  total  computational  load  is  roughly  2 n  additions  plus 
nlogn  operations  for  the  sort. 

2.  Fogel  and  Gavish  have  proposed  the  following  estimate  for  r  [9], 

~  1 

T=  - ; - .  (15) 

argmax/  Y^=i  e27"-^ 

The  denominator  arises  from  spectrum  analysis  of  a  delta-train  function  with  unit  impulses 
at  the  occurrence  times  (see  also  Bartlett  [1]).  This  algorithm  has  a  simple  program  state¬ 
ment  and  requires  no  tuning  parameters  other  than  choice  of  /.  In  our  tests  this  algorithm 
performed  well  given  that  the  percentage  of  missing  observations  was  not  too  large  and  given 
some  a  priori  knowledge  of  a  range  of  /  over  which  to  search.  The  algorithm  requires  fine  steps 
in  /  when  searching  due  to  the  convergence  properties  of  the  sum  of  exponentials.  Without 
this  fine  search  it  is  possible  to  miss  /  =  1/t  altogether.  This  is  an  important  issue  because 
the  number  of  frequency  bins  in  the  search  affects  the  computational  load. 

Ignoring  computation  of  the  exponential  function,  directly  computing  the  denominator  of  (15) 
requires  2 nF  complex  multiplies,  where  F  is  the  number  of  frequency  bins.  Given  that  a  fine 
frequency  sampling  is  required,  the  computational  load  can  grow  significantly  without  some  a 
priori  knowledge  of  where  to  search.  It  is  possible  to  compute  the  denominator  using  an  FFT. 
However,  the  input  to  the  FFT  is  a  delta-train  (of  mostly  zeros)  with  length  determined  by 
the  time  span  of  the  set  S  and  the  system  clock  rate.  Thus,  the  FFT  is  not  of  length  n,  but 
rather  some  significantly  larger  value. 

Use  of  (15)  can  be  contrasted  with  the  proposed  modified  Euclidean  algorithms  that  are 
independent  of  choice  of  r.  As  shown  above  in  Remark  1,  the  robust  single  iteration  algorithm 
has  a  relatively  small  computational  load.  In  addition,  our  algorithms  can  yield  estimates  of 
r  despite  a  very  high  percentage  of  missing  observations.  These  remarks  suggest  that  one 
possible  approach  to  estimation  of  r  is  to  use  the  modified  Euclidean  algorithm  to  obtain  an 
estimate  and  then  fine  tune  this  estimate  using  another  algorithm,  e.g.,  (15). 
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5  Simulation  Results  and  Discussion 


5.1  Simulation  Results 

In  this  section  we  present  simulation  results  and  discuss  performance  of  the  proposed  algorithms. 
All  estimates  and  their  standard  deviations  are  based  on  averaging  over  100  Monte-Carlo  runs. 
The  number  of  data  samples  is  given  by  n  (see  equation  (1)).  Estimates  of  r  are  labeled  t,  with 
experimental  standard  deviations  in  parentheses.  Without  loss  of  generality,  we  take  r  =  1  in 
all  experiments.  This  choice  is  arbitrary,  and  any  real  positive  number  will  yield  similar  results. 
Choice  of  initial  phase  cf)  is  also  arbitrary  as  the  algorithms  are  independent  of  its  choice  due  to 
the  ordering  and  differencing  operations.  The  noise  values  r)j  are  modeled  as  uniformly  distributed 
with  pdf  fv(r])  U[—  y,  y].  So,  for  example,  A  =  10  1  implies  random  phase  jitter  that  is  ±5% 
of  the  period  r  =  1. 

(1)  Noise-free  estimation. 

In  this  example  we  examine  the  effects  of  changing  n  and  the  percentage  of  missing  observations 
on  the  modified  Euclidean  algorithm  of  Section  2.  The  data  are  noise-free,  i.e.,  r)j  =  0  for  all  j. 
In  this  case  the  algorithm  converges  to  the  exact  value  of  r  =  1  with  standard  deviation  equal  to 
zero,  or  in  some  cases  for  small  n  to  some  multiple  of  r.  The  jumps  in  the  kj’s  were  modeled  as 
uniformly  distributed  on  the  (discrete)  interval  [1  ,  M].  Results  are  shown  in  Table  2  where  %miss 
denotes  the  experimentally  determined  average  percentage  of  missing  observations,  and  iter  is  the 
average  number  of  iterations  required  to  converge. 

The  top  half  of  Table  2  illustrates  the  effect  of  changing  M,  and  therefore  changing  the  percent¬ 
age  of  missing  observations.  Given  insufficient  data  the  algorithm  will  converge  to  a  multiple  of  r. 
Columns  labeled  r,  2 r,  etc.,  indicate  the  percentage  of  runs  that  converged  to  these  values.  The 
algorithm  is  able  to  choose  r  correctly  based  on  n  =  10  data  samples,  even  with  99.998%  of  the 
possible  observations  missing.  Convergence  in  the  noise-free  case  depends  on  n  but  is  independent 
of  M,  as  implied  by  the  analysis  of  Section  3.  The  bottom  half  of  Table  2  illustrates  the  effect  of 
changing  n  for  M  fixed.  Reliable  results  are  achieved  for  n  >  10. 

(2)  Noisy  data  and  the  multiple-iteration  algorithm. 
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In  this  example  we  implement  the  modified  Euclidean  algorithm  (with  averaging)  of  Section  4. 
Results  are  shown  in  Table  3.  The  missing  observations  were  modeled  using  a  Bernoulli  process 
with  parameter  A,  A  is  the  noise  parameter,  and  iter  is  the  mean  number  of  iterations  required 
to  converge.  Note  that  A  corresponds  to  the  expected  percentage  of  missing  observations  (see  (14) 
and  discussion). 

The  top  half  of  Table  3  shows  that  estimates  of  r  degrade  as  the  noise  increases,  as  we  would 
reasonably  expect.  The  bottom  half  of  Table  3  shows  the  effects  of  increasing  A,  hence  increasing  the 
percentage  of  missing  observations.  The  performance  is  essentially  unchanged  with  0.6  <  A  <  0.9. 
The  last  entry  in  Table  3  shows  accurate  estimation  of  r  from  100  data  samples  with  ±5%  phase 
jitter  and  90%  of  the  observations  missing.  It  is  possible  to  select  A  >  0.9.  However,  estimation  of 
f  becomes  less  reliable.  For  example,  with  A  =  0.95,  3  out  of  100  trials  resulted  in  poor  estimates 
of  f  that  in  turn  resulted  in  poor  estimates  of  r,  while  the  other  97  trials  produced  estimates  close 
to  r  =  1. 

(3)  Noisy  data  and  the  single  iteration  algorithm. 

Next  we  repeat  Example  2  but  now  employ  the  single  iteration  algorithm  of  Section  4.  Results 
are  shown  in  Table  4.  Here,  stepwidth  is  the  number  of  data  samples  isolated  in  the  lowest  step 
and  whose  mean  resulted  in  f .  The  top  half  of  the  table  shows  the  effect  of  increasing  the  noise. 
As  we  expect  the  estimate  f  degrades  for  higher  noise.  Note  that  the  value  of  stepwidth  reveals 
that  the  algorithm  is  able  to  capture  the  approximately  50  data  samples  expected  in  the  lowest 
step  (due  to  A  =  0.5). 

The  bottom  half  of  Table  4  depicts  the  effects  of  increasing  A.  Here,  estimates  of  r  degrade 
slightly  as  A  increases.  The  degradation  is  due  to  the  decreasing  value  of  stepwidth.  As  A  increases 
toward  1.0  fewer  data  samples  are  available  in  the  lowest  step,  hence  the  mean  of  these  samples 
becomes  a  less  accurate  estimate  of  r.  This  effect  can  be  counterbalanced  by  increasing  n  since  the 
value  of  stepwidth  is  approximately  (1  —  A )n  data  samples. 

(4)  Noisy  data  with  outliers. 

Next  we  test  the  robust  single  iteration  algorithm  of  Section  4  by  adding  outliers,  with  results 
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shown  in  Table  5.  Here,  %outliers  indicates  the  percentage  of  the  n  =  100  data  samples  that 
are  outliers.  These  were  generated  by  choosing  them  to  be  uniformly  distributed  over  the  span  of 
the  set  S.  So,  e.g.,  with  n  =  100  and  %outliers  =  .15,  the  data  set  was  generated  with  85  data 
samples  and  then  15  outliers  were  inserted  for  each  Monte  Carlo  run,  maintaining  the  total  number 
of  samples  at  n  =  100. 

Also  shown  in  Table  5  is  the  step-width  threshold  xo ■  The  last  column,  labeled  %runs ,  indicates 
the  number  of  Monte  Carlo  runs  out  of  100  that  yielded  estimates  of  r.  Failure  of  the  algorithm  to 
produce  an  estimate  occurs  when  there  is  no  step  that  exceeds  the  step-width  threshold  xo-  Note 
the  last  entry  in  Table  5  indicating  successful  estimation  of  t  from  n  =  100  data  samples  with  80% 
missing  observations  and  10  of  the  100  data  samples  being  arbitrary  outliers. 

The  results  of  Table  5  are  somewhat  distorted  by  our  method  of  generating  outliers.  For  small 
A  (fewer  missing  observations)  the  total  time  span  of  the  original  data  set  S  is  considerably  smaller 
than  for  larger  values  of  A.  Thus,  for  small  A,  the  outliers  are  more  closely  spaced  in  time  in 
the  original  data  set  S,  and  therefore  more  likely  to  result  in  a  false  step  with  mean  less  than 
r.  Counterbalancing  this  effect  requires  a  larger  step-width  threshold  xo ■  The  values  for  xq  were 
chosen  to  avoid  selecting  such  false  steps  with  mean  less  than  the  true  value  of  r  =  1.  Hence  the 
relatively  smaller  value  of  %runs  for  small  A  as  shown  in  the  Table. 

5.2  Discussion 

Testing  of  the  modified  Euclidean  algorithm  of  Section  2  in  the  noise-free  case  (Example  1  above) 
reveals  excellent  performance  with  very  few  data  samples,  as  theoretically  predicted  in  Section 
3.  This  performance  is  independent  of  the  percentage  of  missing  observations  but  the  repeated 
subtraction  can  lead  to  very  many  iterations  for  extreme  cases. 

Comparison  of  the  results  of  Examples  2  and  3  reveals  that,  in  general,  the  single  iteration 
algorithm  outperforms  the  multi-iteration  algorithm  by  a  slight  margin.  This  somewhat  counter¬ 
intuitive  result  occurs  due  to  the  averaging  over  each  step  at  each  iteration.  Some  step-widths 
will  be  small,  such  that  the  resulting  average  over  a  particular  step  may  have  a  high  variance  due 
to  the  lack  of  data  samples.  Whereas,  in  the  single  iteration  algorithm  only  a  single  average  is 
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computed  over  the  lowest  step.  Our  results  indicate  that,  if  the  pattern  of  missing  observations 
permits,  the  single  iteration  algorithm  is  likely  to  be  preferable  to  the  multi-iteration  version.  The 
single  iteration  version  requires  less  computation  and  requires  only  two  parameters,  the  gradient 
threshold  go  and  the  step-height  threshold  yo-  Our  experiments  indicate  very  good  robustness  to 
choice  of  go .  Choice  of  yo  depends  on  the  coarse  estimate  of  r  we  called  f ,  which  depends  on  getting 
sufficient  data  in  the  lowest  step  after  the  difference  and  sort  operations,  that  in  turn  depends  on 
the  distribution  of  the  missing  observations  in  the  original  data  set  S. 

The  single  iteration  algorithm  can  also  be  made  robust  with  the  addition  of  one  more  parameter, 
the  step-width  threshold  xq.  Proper  choice  of  xq  requires  some  knowledge  of  the  distribution  of 
the  missing  observations.  However,  proper  choice  can  insure  very  robust  performance  with  heavily 
contaminated  data.  Obtaining  robust  behavior  with  respect  to  outliers  requires  choices  that  are 
not  likely  to  be  optimal  in  all  scenarios.  Our  experimental  choices  are  based  on  experience  and  the 
scenarios  considered. 

Finally  we  note  experimentation  with  using  the  median  rather  than  the  mean  when  averaging 
over  steps,  in  an  attempt  to  further  combat  outliers.  While  the  median  did  produce  slightly  better 
results  in  some  cases,  the  level  of  improvement  was  very  slight. 

6  Final  Remarks 

We  have  developed  modifications  of  the  Euclidean  algorithm  that  are  based  on  subtraction  rather 
than  division,  and  applied  these  algorithms  to  determine  the  period  from  a  sparse  set  of  noisy 
measurements.  These  results  are  based  on  the  model  given  by  (1),  assuming  a  single  period  r. 
The  algorithms  are  computationally  straightforward  and  require  no  a  priori  information.  In  the 
noise-free  case  our  algorithm  is  equivalent  to  the  Euclidean  algorithm.  It  converges  with  very  high 
probability  given  n  =  10  or  more  data  samples.  Modifications  of  the  original  algorithm  work  well 
even  in  the  presence  of  noise  and  outliers. 

It  is  natural  to  consider  an  extension  to  the  case  of  multiple  r’s,  say  in  a  set  {tj}”=1.  We 
ask  what  number  t  will  be  produced  by  the  algorithm.  This  divides  into  two  cases.  The  first 
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case  considers  {Tj}  such  that  the  Tj' s  are  all  rational  multiples  of  each  other.  The  algorithm  will 
produce  a  number  r  such  that  for  each  Tj,  there  exists  a  positive  integer  c3  such  that  CjT  =  Tj. 
The  number  r  is  just  gcd(ri, . . . ,  rn).  (This  more  general  notion  of  gcd  is  consistent  with  Euclid’s 
original  algorithm,  e.g.,  see  Knuth  [19],  pp.  317-320.  Also  see  Proposition  2.1.)  For  example,  if 
Ti  =  5,  T2  =  §,  and  T3  =  |,  then  r  =  This  result  is  easily  seen,  for  it  is  just  as  if  we  applied  the 
algorithm  directly  to  the  set  {Tj}. 

The  second  case  considers  { Tj }  such  that  at  least  two  of  the  Tj' s  are  not  rational  multiples  of 
each  other,  for  example,  t\  =  72  =  \/2.  The  values  of  the  Tj' s  are  incommensurable,  that  is, 

there  is  no  fundamental  unit  in  which  all  of  the  Tj' s  can  be  expressed  as  multiples.  This  follows 
from  the  fact  that  given  a  fixed  irrational  7  and  any  tj  >  0,  there  exists  positive  integers  m  and 
n  such  that  | m  —  rr/l  <  r).  (Also  see  Weyl’s  Equidistribution  Theorem,  Dym  and  McKean  [7],  pp. 
54-56.)  In  this  case,  the  algorithm  will  in  theory  produce  r  =  0,  although  this  is  sensitive  to  the 
noise  thresholding  and  truncation  error. 

Finally,  we  discuss  a  least  squares  procedure  for  fine  tuning  our  estimate  of  r.  Given  the  original 
data  set  S,  estimate  r  and  <f>.  Given  f  and  </>,  get  the  set  {kj}  via  equation  (3)  (see  the  remark  at 
the  end  of  Section  2).  This  last  step  is  sensitive  to  the  quality  of  the  estimates  f  and  4>,  and  to  the 
distribution  of  the  rjj's.  Having  estimated  the  set  {kj},  we  then  form  the  function 

f(t )  =  hT)2i  (!6) 

j 

for  t  >  0,  and  solve  for  r  such  that 

—F(r)  =  -2^kj{Sj-4>-kjT)=Q.  (17) 

i 

Since  d2/dT2F(T)  =  2  J^j(kj)2  >  0)  this  value  of  r  will  minimize  F(t). 
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7  Appendix:  Proofs 


Proof  of  Proposition  2.2  : 

Let  a  be  a  positive  integer  such  that  a\kj  for  j  =  l,...,n.  Then,  a\(kj—kj+\)  for  j  =  1, 

and  Q!|A:n. 

Conversely,  assume  /?  is  a  positive  integer  such  that  (3\(kj  —  kj+ 1)  for  j  =  1, . . . ,  n  —  1,  and  [3\kn. 
Therefore,  there  exists  positive  integers  c  and  d  such  that  c/3  =  kn  and  d(3  =  (&n_i  —  kn).  Thus, 
df3  +  kn  —  (d  +  c)/3  =  kn-i,  and  so  /3\kn-i-  By  complete  induction,  {3\kj  for  j  =  1, . . . ,  n. 

Therefore,  since  the  sets  {kj}  and  {(kj  —  kj+ 1)}  U  {kn}  have  the  same  divisors,  their  gcd’s  are 
equal.  □ 

Remark  :  A  proof  of  gcd(A:i ,  k?)  =  gcd((A;i  —  /C2),  ^2)  with  an  outline  of  the  algorithm. 

We  may  represent  the  Euclidean  algorithm  applied  to  ki,k2,  for  k\  >  k2,  as  follows. 

ki  =  k2  •  qi  +  n  :  0  <  n  <  k^ 
k2  =  r\  ■  q2  +  r2  :  0  <  r2  <  ri 

rk-2  =  rk- 1  -qk  +  rk  ■  0  <rk<  rk-X 

rk- 1  =rk-qk.  (18) 

The  procedure  terminates  when  rk+x  —  0.  This  gives  gcd(ki,k2)  =  rk.  Since  qk  >  1, 

ki-k2  =  k2-  ( q\  -  1)  +  n  .  (19) 

We  may  replace  the  top  line  of  (18)  with  (19),  yielding  the  same  result.  Thus, 

gcd(Aq,  A:2)  =  gcd((fci  -  k2),k2).  (20) 

Therefore,  given  two  positive  integers  k\  and  k2,  kx  >  k2,  we  can  reformulate  the  Euclidean 
algorithm  using  subtraction  rather  than  division.  This  is  stated  in  pseudocode  as  follows. 
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FUNCTION  gcd(fci,ib2) 

Begin 

Do  while  &2  ^  0 
k\=k\-  fc2 

If  h<k2,  SWAP(fci,  A:2) 
Enddo 

gcd(fci,  &2)  =  ki 
End 


Proof  of  Proposition  2.3  : 

Let  a  be  a  positive  integer  such  that  a\(kj  —  kj+ 1)  for  j  =  1, . . . , n  —  1.  Therefore,  a\(kn-i  —  kn) 
and  a\(kn-2  —  kn-i),  and  so  there  exist  mi,  m2  such  that  kn- \  —  mia  +  kn  and  kn-2  =  m2a  +  /cn_i. 
Substituting,  we  get  that  /c„_2  =  (m2  +  mi)a  +  kn ,  and  so  a\(kn-2  —  kn).  By  complete  induction, 
we  get  that  a\(kj  —  kn),  for  j  =  1, . . .  ,n  —  1. 

Conversely,  assume  (5  is  a  positive  integer  such  that  (3\(kj  —  kn)  for  j  =  1, . . .  ,n  —  1.  Thus 
0\{kj  —  kn)  and  fi\(k]+i  —  kn),  and  so  there  exists  mj,  mry\  such  that  kj  —  rrijfl  +  kn  and  k^+i  = 
m,j+ \f3  +  kn.  Thus,  kj  -  kj+i  =  ( rrij  —  rrij+ 1)/3,  and  so  (3\(kj  -  kj+ 1). 

Therefore,  since  the  sets  {(kj  —  kj+ 1)}  and  {(kj  -  kn)}  have  the  same  divisors,  their  gcd’s  are 
equal.  □ 


Proof  of  Theorem  3.2  : 

Let  [rcj  denote  the  floor  function  of  x,  namely 


\x\  =  :  k  €  Z}  . 


(21) 


Recall  that  Nn(£)  =  card{(A:i, . . . ,  kn)  E  {l,...,£}n  :  gcd(^i, . . . ,  kn)  =  1}  is  the  number  of 
relatively  prime  elements  in  {1, . . . ,  £}n.  We  claim  that 


E 

Pi<Pj 


L— J 

k  Pi  ■  Pi  , 


-  E 

Pi<Pj<Pk 


L: 


Pi  ■  Pj  ■  Pk 


j  + 


(22) 


This  is  seen  as  follows.  Choose  a  prime  number  pi.  The  number  of  integers  in  {1  such 

that  pi  divides  an  element  of  that  set  is  (Note  that  is  is  possible  to  have  pi  >  £,  because 

[jfj  =  0.)  Therefore,  the  number  of  n-tuples  (ki, . . . ,  kn)  contained  in  the  lattice  {1, . . . ,  £}n  such 
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that  pi  divides  every  integer  in  the  n-tuple  is 

(£)'• 

Next,  if  pi  ■  pj  divides  an  integer  k ,  then  pl\k  and  pj\k.  Therefore,  the  number  of  n-tuples 
kn)  contained  in  the  lattice  {1, . . . ,  £}n  such  that  pi  and  pj  both  divide  every  integer  in  the 
n-tuple  is 


where  the  last  term  is  subtracted  so  that  we  do  not  count  the  same  numbers  twice  (in  a  simple 
application  of  the  inclusion-exclusion  principle). 

Continuing  in  this  fashion,  for  three  integers,  say  pt  <  pj  <  pk,  the  number  of  n-tuples 
(ki, . . . ,  kn)  contained  in  the  lattice  {1, . . . ,  £}n  such  that  pi,  pj,  and  pk  all  divide  every  integer 
in  the  n-tuple  is  given  by  the  inclusion-exclusion  principle  as 


(25) 

We  can  therefore  see  by  induction  that  the  number  of  n-tuples  {h,... , kn)  contained  in  the 
lattice  {l,...,f}n  such  that  Pi,Pj,pk,  ■  ■  ■  all  divide  every  integer  in  the  n-tuple  is  given  by  the 
inclusion-exclusion  principle  as 


L — —  J 

Pi  •  Pj  ■  Pk 


(26) 


But  this  counts  the  complement  of  Nn(£)  in  the  lattice  {1, . . . ,  £}n.  Therefore, 

w«M  =  <n-£(l^j)"+  E  E  (*) 

Pi  'Pi  '  Pi<Pj  \  Pi  Pj  )  pi<Pj<Pk  \  Pi  Pj  Pk  J 


Next,  we  observe  that 
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where  the  last  sum  over  j  G  N  \  {1}.  Since  n  >  2,  this  series  is  convergent.  Therefore,  each  term 
in  the  expansion  of  Nn^  is  convergent. 

We  recall  the  Mobius  function  p,  which  is  defined  as  follows. 


Mi)  =  i, 


Mm)  =  { 


0  if  m  is  divisible  by  the  square  of  a  prime , 

(— l)r  if  m  =  pi  •  P2  •  •  •  •  •  pr,  where  pi,P2,  ■  ■  •  ,Pr  are  all  distinct  primes . 


(28) 


Euler  showed  that 


1  +  P5,  (Pi-Pj)n 


p;- 

Pi 


=  E —?  =  [«")]“’ 


m" 


where  the  last  sum  is  over  m  G  N.  For  n  >  2,  Em  Em  are  absolutely  convergent.  Thus, 
the  equality  Em  =  [C(n)]_1  follows  because  for  j,  k,m,n  G  N, 

E^E^  =  E^  =  EsEMd)  =  i, 


m  mn  .  j 

m  j  J 


kr 

k  d\k 


where  we  have  used  the  fact  that  both  series  in  the  first  term  converge  absolutely  and  thus  can 
be  multiplied  together  by  adding  all  possible  products  of  terms  and  rearranging  in  any  order  (see 
Ireland  and  Rosen  [14]  and  Leveque  [21]). 

Let 


Ms? 

with  sum  over  j  €  N  \  {!}•  Since  n  >  2  and  the  sum  is  over  j£N  \  {1}, 


(29) 


(30) 
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Since  the  A;th  term  in  the  expansion  of  pP-  is  dominated  by  M*,  and  Y,k  Mib  is  convergent,  we  may 
apply  the  Weierstrass  M  test,  and  evaluate  term-by-term.  We  use  the  fact  that  x  —  l  <  \x\  <  x 
for  all  x,  and  so  lim^oo  “  =  1.  Therefore,  we  have 


*->oo  in 


=  lim (V- W|AlV+  E  (l— jV-  E  fi  — 

^  ^\pi  J  p%  -Pp)  Pi£<Pk  VLw  • 


Pj  ■  Pk 


J  +  ... 


-  1  -  E^r  +  E 


E 


Pi 


Pi  P%  (Pi  •  Ppn  Pi<ip<Pk  (Pi  ‘  Pj  ■  Pk) 


+  ... 


y 

mn 

m 

=  [C(n)]-1, 


where  the  last  sum  is  over  me  N.  This  completes  the  proof  of  the  theorem.  □ 

Proof  of  Corollary  3.1  : 

By  Proposition  2.3, 


gcd((Aq  ^2))  (^-2  ^3)1  •  •  •  >  (^n— 1  ^n))  — gcd((Aq  kn),  (^2  ^n)i  ■  •  •  i  {kn— 1  &n))  • 


Thus,  computing  the  gcd  of  the  left  hand  side  is  the  same  as  computing  the  gcd  of  the  set 
{ki,k%, . . . ,  kn-\}  shifted  by  kn.  Now,  since  kj  —  kj+\  >  1,  we  see  that  we  can  repeat  the  analysis 
of  Theorem  3.2  in  the  (n  —  1) -dimensional  sublattice  of  N"  consisting  of  the  of  (n  —  l)-tuples  of 
the  form 

{(m  +  kn), . . . ,  (m  +  kn)},  1  <  m  <  l . 


In  this  sublattice,  for  sufficiently  large  t,  the  number  of  elements  that  are  divisible  by  some  given 
prime  p  is  essentially  equal  to  the  number  of  elements  in  {1, . . .  , divisible  by  p.  Moreover, 
if  we  denote  the  number  of  relatively  prime  (n  -  l)-tuples  in  this  shifted  lattice  by  Nn-i(£a),  we 
have  that 


lim 

*— >oo 


Nn-i(is) 

£n-l 


lim 

*— >oo 


iVn-l(l) 

fn-\ 


(Asymptotically,  the  fact  that  the  sublattice  is  shifted  by  some  finite  fixed  amount  will  not  matter.) 
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But  then,  by  Theorem  3.2, 


lim 

i—yoo 


Nn-i(i) 

(n-\ 


^[an-1)]-1 


•  □ 


Proof  of  Proposition  3.1  : 

Since 

00 

CM  =  Y71^ 

71  =  1 

and  lj  >  1, 


i  <  CM  =  i  +  —  +  —  +  —  +  —  +  ... 

.  ,  1  1  1  ii  i 

-  1  +  ^  +  2^  +  ^  +  ---  +  ^  +  gJ  +  ---  +  8^  +  --- 

s - - - '  ' - v - ' 

4-times  8— times 

~  1  -  £  _  1 -21-"  ' 

As  u)  — >  oo,  (1  —  21_w)  — >  1+.  Therefore,  by  the  Squeezing  Theorem, 

[CM]-1 — >l~asw — >oo.  D 
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Figure  1 


j 


1.  Plot  of  example  data  set  after  one  iteration  of  the  modified  Euclidean  algorithm  of  Section  2. 
The  data  is  sorted  in  descending  order  into  steps  centered  around  multiples  of  r  (r  =  1  in  this 
example).  The  stepwidths  are  a  function  of  the  distribution  of  the  kj’s  in  the  original  data  set  S. 
The  lowest  (rightmost)  step  is  centered  around  the  true  value  of  r  =  1. 
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Table  1:  Some  values  of  the  Zeta  function  ((n)  and  l/((n). 


n 

C(«) 

i/C(») 

2 

1.64493 

0.6079 

4 

1.08232 

0.9239 

6 

1.01734 

0.9830 

8 

1.00407 

0.9959 

10 

1.00099 

0.9990 

12 

1.00024 

0.9998 

16 

1.00001 

1.0000 

Table  2:  Results  from  simulation  Example  1,  noise-free  estimation  of  r  with  the  modified  Euclidean 
algorithm. 


n 

M 

%miss 

iter 

T 

2  T 

3  r 

>  3  T 

10 

101 

81.69 

3.3 

100% 

0 

0 

0 

10 

102 

97.92 

10.5 

100 

0 

0 

0 

10 

103 

99.80 

46.5 

100 

0 

0 

0 

10 

104 

99.98 

316.2 

100 

0 

0 

0 

10 

105 

99.998 

2638.7 

100 

0 

0 

0 

4 

102 

97.84 

15.2 

82% 

12 

4 

2 

6 

102 

97.81 

14.2 

97 

3 

0 

0 

8 

102 

97.96 

10.2 

98 

1 

1 

0 

10 

102 

97.95 

10.2 

99 

1 

0 

0 

12 

102 

97.95 

8.6 

100 

0 

0 

0 

14 

102 

97.97 

7.4 

100 

0 

0 

0 

Table  3:  Example  2  results,  estimation  of  r  from  noisy  measurements  using  the  modified  Euclidean 
algorithm  (with  averaging). 


n 

A 

A 

t  (std) 

iter  (std) 

100 

.5 

10"3 

1.0000  (.0001) 

3.5  (.6) 

100 

.5 

10-2 

1.0000  (.0014) 

3.5  (.5) 

100 

.5 

10-1 

1.0002  (.0169) 

3.5  (.6) 

100 

.5 

2  x  10_1 

1.0032  (.0212) 

3.5  (.5) 

100 

.6 

10”1 

0.9998  (.0158) 

3.7  (.6) 

100 

.7 

10"1 

0.9997  (.0200) 

4.0  (.6) 

100 

.8 

io-1 

1.0002  (.0202) 

3.9  (.4) 

100 

.9 

io-1 

1.0038  (.0213) 

4.2  (.6) 
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Table  4:  Example  3  results,  estimation  of  r  from  noisy  measurements  using  the  single  iteration 
algorithm. 


n 

A 

A 

f  (std) 

stepwidth  (std) 

100 

.5 

10-3 

1.0000  (.00001) 

49.0  (4.6) 

100 

.5 

10~2 

1.0000  (.0004) 

48.3  (4.5) 

100 

.5 

10-1 

0.9994  (.0043) 

49.2  (5.0) 

100 

.5 

2  x  lO-1 

0.9959  (.0104) 

47.1  (5.2) 

100 

.6 

HT1 

1.0000  (.0056) 

38.5  (4.9) 

100 

.7 

ur1 

1.0009  (.0061) 

28.3  (5.0) 

100 

.8 

HT1 

0.9986  (.0089) 

18.6  (3.9) 

100 

.9 

io-1 

0.9988  (.0141) 

8.1  (2.7) 

Table  5:  Example  4  results,  for  noisy  sparse  data  with  outliers,  using  the  robust  single  iteration 
algorithm. 


n 

A 

A 

%outliers 

f  (std) 

x0 

%runs 

100 

.05 

IO-1 

0.15 

0.9983  (.0035) 

35 

93 

100 

.05 

io-1 

0.25 

0.9973  (.0054) 

35 

86 

100 

.05 

IO"1 

0.35 

0.9928  (.0099) 

35 

52 

100 

.10 

io-1 

0.05 

0.9994  (.0027) 

25 

100 

100 

.10 

10~x 

0.15 

0.9991  (.0036) 

25 

99 

100 

.10 

io-1 

0.25 

0.9975  (.0062) 

25 

98 

100 

.10 

io-1 

0.35 

0.9958  (.0105) 

25 

90 

100 

.10 

io-1 

0.45 

0.9839  (.0149) 

20 

16 

100 

.25 

IO'1 

0.05 

1.0006  (.0032) 

15 

100 

100 

.25 

io-1 

0.15 

0.9987  (.0037) 

15 

100 

100 

.25 

10"1 

0.25 

0.9987  (.0052) 

15 

100 

100 

.25 

io-1 

0.35 

0.9953  (.0095) 

15 

93 

100 

.50 

io-1 

0.05 

1.0014  (.0048) 

10 

100 

100 

.50 

io-1 

0.15 

1.0004  (.0069) 

10 

100 

100 

.50 

io-1 

0.25 

1.0003  (.0071) 

15 

100 

100 

.50 

io-1 

0.35 

1.0004  (.0110) 

15 

98 

100 

.75 

IO"1 

0.05 

1.0043  (.0092) 

5 

100 

100 

.75 

HT1 

0.10 

1.0036  (.0344) 

10 

100 

100 

.75 

io-1 

0.15 

1.0064  (.0189) 

10 

100 

100 

.75 

HT1 

0.20 

1.0004  (.0211) 

10 

100 

100 

.80 

io-1 

0.05 

1.0038  (.0103) 

5 

100 

100 

.80 

IO-1 

0.10 

1.0061  (.0143) 

5 

100 
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